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ABSTRACT 

The intracluster medium (ICM) has been suggested to be buoyantly unstable in the presence of 
magnetic field and anisotropic thermal conduction. We perform first cosmological simulations of 
galaxy cluster formation that simultaneously include magnetic fields, radiative cooling and anisotropic 
thermal conduction. In isolated and idealized cluster models, the magnetothermal instability (MTI) 
tends to reorient the magnetic fields radially whenever the temperature gradient points in the direction 
opposite to gravitational acceleration. Using cosmological simulations of the Santa Barbara cluster 
we detect radial bias in the velocity and magnetic fields. Such radial bias is consistent with either 
the inhomogeneous radial gas flows due to substructures or residual MTI-driven field rearangements 
that are expected even in the presence of turbulence. Although disentangling the two scenarios is 
challenging, we do not detect excess bias in the runs that include anisotropic thermal conduction. 
The anisotropy effect is potentially detectable via radio polarization measurements with LOFAR 
and Square Kilometer Array and future X-ray spectroscopic studies with the International X-ray 
Observatory. We demonstrate that radiative cooling boosts the amplification of the magnetic field 
by about two orders of magnitude beyond what is expected in the non-radiative cases. This effect 
is caused by the compression of the gas and frozen-in magnetic field as it accumulates in the cluster 
center. At z = the field is amplified by a factor of about 10^ compared to the uniform magnetic 
field evolved due to the universal expansion alone. Interestingly, the runs that include both radiative 
cooling and anisotropic thermal conduction exhibit stronger magnetic field amplification than purely 
radiative runs, especially at the off-center locations. In these runs, shallow temperature gradients 
away from the cluster center make the ICM neutrally buoyant. Thus, the ICM is more easily mixed in 
these regions and the winding up of the frozen-in magnetic field is more efficient resulting in stronger 
magnetic field amplification. We also demonstrate that thermal conduction partially reduces the gas 
accretion driven by overcooling despite the fact that the effective conductivity is suppressed below the 
Spitzer-Braginskii value. 

Subject headings: conduction - cooling fiows - galaxies: clusters: general - galaxies: active - instabil- 
ities - X-rays: galaxies: clusters 



1. INTRODUCTION 

Thermal conduction may play an important role 
in the evoluti on of the intracluster medium (e.g., 
iVoigt fc FabianI ()2664l )). Sufficiently strong conduction 
can offset radiative cooling losses in massive clusters 
and reduce the energy requirements on active galac- 
tic nuclei (AGN) feedback that is required to prevent 
overcooling in less massiv e clu sters and grou ps (see 
iMcNamara fc NulsenI (pOOl and E orma 3 (l20Tnl) for re- 
views). It has also been suggested that not only can 
thermal conduction serve as a mechanism for cool core 
heating but that it i s very important for the stabi l- 
ity of the s e syst ems (jRuszkowski fc BegelmanI ()2002D . 
iGuo et all pOOSi ). Ruszkowski & Oh, in prep.). It may 
also be responsible for setting a critical central entropy 
threshold belo w which star formation is possible in clus- 
ter cool cores (jVoit et al.l 120081 ). The recently discov- 
ered bim odality in the distribution of the cluster centra l 
entropy (|Cavagnolo et al.l 120091: [Sanderson et al.l |2009() 



may be due to the combination of AGN feedback from 
the brightest cluster galaxies that stabilizes low entropy 
clusters and a combination of mergers and thermal con- 
duction th at stabilize higher centr al en t ropy clusters 
(IGuo et al. 2008: Ruszkowski & Oh 2 010bl : lParrish et al.l 
|2010( ). As" such, thermal conduction may be important 
for understanding of the feeding of the most massive 
black holes in the Universe. In cluster outskirts ther- 
mal c onduction may flat ten the temperature distribu- 
tions (jParrish et al.H2008D . which may have consequences 
for the cluster mass estimates. This may have possible 
impact for precision cosmology as it relies on accurate 
mass measurements in the most massive clusters. 

As a plasma transport process, thermal conduction 
is closely linked to gas viscosity. Both types of 
transport processes may explain various X-ray obser- 
vatioi is. For exaniple, r ecent Chandra observations of 
M87 (jWerner et al.ll2010l ) show that the temperature in 
the shells centered on the cluster center is remarkably 
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isothermal. They suggest that such a high degree of 
isothermaUty is consistent with effective heat conduction 
in the tangential direction. Moreover, they also attribute 
the presence of small scale metallicity gradients to rela- 
tively weak level of small scale turbulence, which could 
be consistent with viscous damping of gas motions. AGN 
are known to generate intermitt ent outflows c ausin g 
weak shocks a n d sou n d waves (e.g.. iFabian et al.l ()2003[ ). 
iForman et all (|2005D . iFinoguenov et al.l (|2008[) V The 
dissipation of the energy contained in these waves may be 
sufficient to offset radiative cooling of the gas. It has bee n 
shown obs ervationally (in Perseus: Fab ian et al.l (|2003f ): 
in Virgo: iForman et al.l ([20()5l)) and u sing numerical 
simula tions ([Ruszkowski et al.l l2004al lbl: iBriiggen et al.l 
l2005aD that Spitzer-Braginskii viscosity and/or conduc- 
tion is sufficient to dissipate such waves and heat the 
ICM efficiently. Moreover, the compact morphology of 
the buoyantly rising AG N is consistent with the pres- 
ence of effective viscosity (|Revnolds et al.li2005 ,') although 
it co uld also be explained by the magnetic draping ef- 
fect dR obinson et al."2004!: ' Ruszkowsk i et al.!l2007ll2008l: 
IDursi fc Pfrommcr 2008; O 'Neih et al. 2009). Other- 
wise, the bubbles would be disrupted by Ray leigh- Taylor 
instabilities. Gas viscosity has also been shown to play a 
role in shaping the pr operties of cold fronts caused b y the 
sloshing motion (e.g.. iZuHone fc MarkevitchI (|2009[ )V 

The simple picture of viscosity and conduction is com- 
plicated by the presence of magnetic fields that are 
known to be present in the ICM (Enfilin & Vogt 2003,; 
IVogt fc EnfiimI l200l iFeretti fc Johnston-HoUitti 120041 ) . 
Magnetic fields suppress thermal conduction in the di- 
rection perpendicular to the B-field. However, even in 
the case of highly tangled magnetic fields, the effec- 
tive thermal conduction can be a substantial fractio n 
the Spitzer conductivity (jNaravan fc MedvedevI |2"001D . 
The anisotropy due to magnetic fields leads to new phe- 
nomena with potentially important consequences for the 
ICM It has been demon strated both analyticall y (jBalbusI 
l2000t) and numerically (jParrish fc Ston(j|2005l ) that the 
ICM is unstable in the presence of a weak magnetic field 
and anisotropic thermal conduction when the temper- 
ature is increasing in the direction of gravity (magne- 
tothermal instability, MTI). When magnetic fields are 
partially aligned with the temperature gradient, the heat 
flow from hotter to cooler regions along the field lines 
makes such regions more buoyant. This causes the mag- 
netic field to become more aligned with the temperature 
gradient and leads to the instability and preferentially 
radial magnetic fields. A similar instability occurs when 
anisot ropic transport of energy via cosmic rays takes 
place ()Chandran fc DennisI I2006D . More recent analy- 
sis shows that the gas is also unstable when the tem- 
perature decreases in the direction of gravity i n the pres- 
ence of a background heat flux (HBI instability; 'Quataer"^ 
(|2008i) ) which has also been verifi ed by numerical simu- 
lations ()Parrish fc QuataertI 120081 ). The saturated state 
of the HBI corresponds to the magnetic fields oriented 
in the direction perpendicular to gravity. Such field con- 
figuration implies effectively vanishing thermal conduc- 
tion from the hotter outer cluster layers to their cool 
cores. This accelerates the effective cooling rate in the 
core. However, the exact topology of the magnetic fields 
depends also on whether externally imposed turbulence 
driving is present. Such turbulence may come from AGN 



outbursts, galaxy motions and structure formation (ma- 
jor and mi nor mergers). Th is has b een r ecently inves- 
tigate d by iRuszkowski fc OhI ()2010b[) and iParrish et all 
(|2010( ) who showed that there exists a critical level of 
turbulence above which the field can be randomized and 
the conductive heating to the core restored. Although 
this analysis was performed for the HBI, similar argu- 
ments apply to the MTI. 

The exact level of turbulence in the ICM is not 
known but is expected to vary throughout the evolu- 
tion of a cluster . Ho wever, indirect measurements (e.g., 
iChurazov et~an (j2003. )) and upper limits (Sander s et al.l 
(I2010D ) are available. For example, some observa- 
tional estimates put the ICM velocity at a level of lo- 
cal so und speed (jMarkevitch et al.l 120021 : iMahdavi et al.l 
120071). while others suggests relatively "calm" ICM 
(jWerner et al.l l2010). Thes e levels o f turb ul ence are ex- 
pected theore tically (e.g. , lEvrardl (|1990f ). iNagai et al.l 
(2003) ■ I Vazza et al.l (|2009[ )'). Future measurements with 
the International X-ray Observatory (IXO) will help to 
determ ine the level of turbule nce in clusters more pre- 
cisely (jHeinz fc BriiggeiJ l2009f ) . Therefore, ideally, we 
need to resort to ab initio cosmological simulations to 
include the effects of structure formation. Using this ap- 
proach we can not only compute the level of the effective 
thermal conduction in the presence of magnetic fields but 
we can also simulate the growth of magnetic field. Both 
MTI and HBI are expected to amplify the fields due to 
kinematic dynamo action although the efficiency of this 
process is modest. However, trapping of gravity modes 
may lead to vorticity gro wth (|Lufkin et al.lll995.) a nd fur- 
ther field amplification ([Ruszkowski fc Ohll2010al ). 

Here we present first cosmological simulations of 
cluster formation that simultaneously include radia- 
tive cooling, magnetic field and anisotropic thermal 
conduction. These simulations are a natural ex- 
tension of our previous work on the role of ra- 
diative cooling, magnetic fields, conduction and vis- 
cosity with the FLASH and ATH ENA codes that 
focused on isolated cool cores (iRuszhowski et alj 
2004arbl: IBriiggen et al.l I2q05at iB ogdanovic et al. 200 J 
Ruszkowski fc OhI l2010bl lal: IParrish et al.l I20ia) . This 
work also builds on previous theoretical and nu- 
merical efforts to simulate the growth of magnetic 
field in cluster formation simulations by our group 
(jBriiggen et all 2005b), with t he F LASH code and 
other teams ('Dolag.eLaD 11999', '2002") with GADGET; 
(CoUins et al. 2009) Enzo^il^i et al. 2008) GosmoMHD, 
(Dubois fc Tcvssicr 2008), RAMSES as weU as the work 
bv .Siiacki fc Springel (20061) . who consid ered non-MHD 
simula tions but included visc osity, and lJubelgas et al.l 
(|2004[ ) and iDolag et all (j2004[ ) who considered conduc- 
tion in hydrodynamical SPH simulations. We also note 
that the field amplificati on may result from purely kinetic 
plasma processes (e.g., ISchekochihin fc CowlevI p007.) . 



iSchekochihin et al.l (|2010[ ). iKunz et al.l (|2010D ) and that 
the large scale turbulence may serve as the energy in- 
put for these mechanisms. However, these processes are 
beyond the scope of this investigation and are not con- 
sidered here. 

The paper is organized as follows. In the next section 
we discuss the methods and the code used for the sim- 
ulations. In Section 3 we discuss our results. Section 
4 presents the conclusions. Appendix A discusses the 
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comparison of the anisotropic conduction module with 
the linear theory predictions. In Appendix B we briefly 
discuss the tests of the implementation of the cosmolog- 
ical terms in the MHD equations. 

2. METHODS 

We solve the following set of MHD equations aug- 
mented by the cosmological expansion terms. 



|+V.(pv)=0 



(1) 



+ V • (pvv - BB) + Vp = -pV(f> - 2-pv (2) 
ot a 



dE 
'dt 



V • (v(S + p) - B(v ■ B)) = -pV0 • V 



+ - [(37 - l)pe + 2pv^] + aC(rph) - V • Fph (3) 



— + V ■ (vB - Bv) = 0, 



where 



P = Pth 



pv 



A^^AnGa-^ip-ip)), 



(4) 

(5) 

(6) 
(7) 



where (p) is the comoving mean density and where pth is 
the gas pressure and e is the gas internal energy per unit 
volume. We assume the adiabatic index 7 = 5/3 and 
the mean molecular weight p = 0.5 in the equation of 
state. In equation (3), C represents the cooling rate per 
unit volume. We use st andard tabulated and pub licly 
available cooling curves (jSutherland fc Dopitalll993l ) for 
metallicity Z = 0.3^0. Our simulations do not include 
star formation or AGN feedback and we set a floor in the 
cooling function at 0.01 keV in physical units. 

The anisotropic thermal conduction heat flux Fph is 
given by 



ph 



-iK(rph)eB(es • VTph), 



(8) 



also imposed an upper limit on conduction such that 
K = a~^K/{cvp), where = (7 - l)M"^prot/fcboitz (the 
extra factor of a comes from the prefactor in Equation 8). 
We set K < if max = 5 X lO^^cm^s"-^. This ceiling was 
introduced in order to prevent extremely small diffusive 
timesteps. We found that changing i^max to larger values 
did not affect our results significantly. Only a relatively 
limited volume far from the cluster center was subject to 
this upper limit. 

We also note that the effective conduction may be 
limited i n the o uter regions of the cluster for physical 
reasons (|Medvede v 2007). Moreover, around or beyond 
the virial radius, the electron-proton energy equilibra- 
tion timescales start to exceed the dynamical or buoy- 
ancy timescale. In this region, the MHD approxima- 
tion starts to break down and the MTI growth rates 
are reduced. However, at smaller radii, the equilibra- 
tion timescale due to Coulomb collisions are comparable 
to or shorter than the buoyancy timescale, some ther- 
mal coupling between electrons and protons begins to be 
possible, and the MHD analysis applies at least approxi- 
mately. Finally, there are independent observational ar- 
guments for fas t equilibration in the ICM. Specifically, 
iMarkevitch fc V ikhlinin (2007,) detected a sharp electron 
temperature jump in the post shock region in the Bullet 
Cluster and used this observation to suggests that the 
electron-proton equilibration is much shorter than the 
collisional timescale. The inferred equilibration timescale 
in the postshock ICM is at least 5 times shorter than the 
Coulomb equilibration time and may even be consistent 
with being instantaneous. While the issue of equilibra- 
tion is clearly an open research topic, here we work under 
the symplifying assumption that equilibration is instan- 
taneous and that the MHD approximation is appropriate, 
and then investigate the consequences of these assump- 
tions. 

The most stringent limitation on the timestep dt in 
the simulation was due to the diffusive process with 
dt ^ {AxY /K. The maximum resolution that we could 
afford computationally was 7 levels of refinement for 
blocks consisting of 16 zones on a side. That is, the 
effective resolution was ~ 31/i~^kpc. The strong limit 
imposed on the timestep can be avoided by employing 
implicit integration methods. We are now implementing 
such methods in the FLASH code (Lee et al., in prep.). 
This approach will significantly speed up the computa- 
tions and will allow us to perform simulations for a range 
of cluster masses (Ruszkowski et al., in prep). 

Equations 1 through 8 were obtained by starting from 
the MHD equations and applying cosmological expansion 
transformation for the spatial gradients and the following 
variable transformations 



where is a unit vector pointing in the direction of 
the magnetic field and k is the Spitzer-Braginskii con- 
duction coefficient given by k{T) = 4.6 x 10~^r^/^erg 
s~^cm~^K^^. In the above equation and in all equa- 
tions below all variables with t he "ph" subscript denote 
physical quantities. Following 'Cowie fc McKed ()1977l ) 
we included the effect of conduction saturation whenever 
the characteristic lengthscale associated with the tem- 
perature gradient exceeds the mean free path, though in 
the bulk of the ICM this effect is not significant. We 



P = apph 
T = a'^Tph 
pe = a/OphCph 



(9) 
(10) 

(11) 
(12) 
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B = fli/^Bpj^, (13) 



where a is the cosmological expansion factor. The 
velocity vector v = x, where x — a~^r is the code 
position vector and r is the physical position vector. 
These transformations together with Equations 1 to 7 
lead to correct scalings of all quantities with the expan- 
sion factor. For example, Tph oc and i3ph oc 
which ensures the conservation of the magnetic flux. 
Note that the variable transformation adopte d here is 
somewhat different than, e.g., that in iLi et all (120081 ) in 
the CosmoMHD code, but the scalings of the physical 
quantities with redshift is correct in both cases. The 
cosmological terms were included using the operator 
splitting technique. Specifically, we computed the up- 
dates to all variables due to the cosmological expansion 
terms by finding the exact solutions to a modified set of 
Equations 1 to 4 that retained only the time derivative 
terms. Tests of the implementation of the cosmological 
terms are presented in Appendix A. 

We used publicly avail able initial condi tions for 
the Santa Barbara cluster (jFrenk et al.l 119991) . These 
conditions correspond to initial Sa density pertur- 
bations spread over 10 Mpc, Hubble constant of 50 
km/s/Mpc and flat geometry in a matter dominated 
universe (A = 0) and result in the formation of a 
cluster characterized by r^-„ ~ 2.7Mpc. This choice 
of publicly available initial conditions allows for stan- 
dardized comparisons with other codes. We employed 
periodic boundary conditions. We note that the initial 
conditions are not sensitive to the initial temperature 
and strength of the magnetic fields as long as they are 
small. The initial magnetic field was assumed to be 
constant in space. While in isolated cluster simulations 
the initial topology of the magnetic field can affect the 
instability growth rate, the frozen-in magnetic field in 
the cosmological simulations is randomized by structure 
formation motions well before the hot virialized and 
relaxed structures are formed. Since early on thermal 
conduction is not expected to be very important, the ini- 
tial field topology is not expected to play a crucial role. 
In fact, it has been demonstrated using cosmological 
MHD simulations without transport processes that the 
results were indeed not sensitive to the exact t opology 
of the magnetic field in the s tatistical sense (Dola g et al.l 
l2002t iBriiggen et al.l l200"5bl) . Therefore, for simplicity 
we assume initial magnetic field that is pointing in 
the same direction and is uniform. We set the initial 
strength of the magnetic field to 10~^^/(47r)^/^ Gauss, 
i.e., the initial physical magnetic field strength was 
10-"/(47r)(l -I- z)i/2 Gauss at the initial redshift of 
z = 20. 

The simulations were performed with the adaptive 
mesh refinement (AMR) FLASH code. FLASH is a 
publicly available code that was in part developed by 
the DOE-supported ASC/ Alliance Center for Astro- 
physical Thermonuclear Flashes at the University of 
Chicago. It is a modular, parallel simulation code 
capable of handling general compressible flow problems 
found in many astrophysical environments. The code 
is parallelized using Message-Passing Interface (MPI) 



library and the IIDF5 or PnetCDF library for parallel 
I/O to achieve portability and scalability on a variety of 
different parallel computers. 

The new d irectionally unsplit st aggered mesh M HD 
solver (USM; ILee fc Deanj ([20091 ): ILee et al.l (pOOflh 'l is 
based on a finite-volume, high-order Godunov method 
combined with a constrained transport (CT) scheme 
which ensures divergence-free magnetic fields. Tests 
of the module demonstrate that it is very robust and 
significantly outperforms the previously implemented 
MHD 8-wave solver in FLASH. 

We implemented the anisotropic conduction uni t 
following the approach of iSharma k, Hammetti ()2007fl . 
More specifically, we applied monotonized central 
(MC) limiter to the conductive fluxes. This method 
ensures that anisotropic conduction does not lead 
to negative temperatures in the presence of steep 
temperature gradients. We verified that the mod- 
ule predicts correct magnetothermal instability growth 
rates. Details of these tests are discussed in Appendix B. 

3. RESULTS 

3.1. Temperature distribution 

Figure 1 shows cross sections through the tempera- 
ture distribution in the cluster. All panels correspond to 
32/i^^Mpc on a side. The minimum and maximum val- 
ues of the temperature are the same in all panels. Top 
row is for the non-conductive cases while the bottom row 
is for the runs with anisotropic thermal conduction. Left 
columns is for non-radiative simulations and the right 
one for the runs with radiative cooling. The upper right 
panel clearly shows a well-developed cool core. Compari- 
son of the top and bottom rows shows that conduction is 
efficient in smearing out the fine structure details in the 
temperature distribution. The cool core in the conduc- 
tive case is partially heated by thermal conduction from 
hotter outer layers of the cluster but it does not become 
isothermal. 

In Figure 2 we show the temperature profiles of the 
cluster for all four cases shown in Figure 1. The profiles 
are not density- weighted. Dashed green curve is for the 
non-radiative run, solid dark blue for anisotropic conduc- 
tion, dashed light blue for radiative cooling, and solid red 
is for radiative cooling with anisotropic conduction. The 
adopted bin size was ~ 31/i~^ kpc. For larger bin sizes, 
the comparison between the temperature profiles and the 
temperature maps shown in Figure 1 is more difficult; 
other one-dimensional profiles presented below use a bin 
size of 125/i~^ kpc, which is especially important for the 
observables inferred from the vector quantities as they 
tend to exhibit larger fiuctuations. The effect of ther- 
mal conduction on the temperature profiles is mild but 
noticeable. This is partially due to the fact that mag- 
netic fields reduce the effective level of conduction be- 
low the Spitzer value . The effect is more pronounced in 
lJubelgas et al.l (|2004[ ) who use full isotropic Spitzer con- 
ductivity in cosmological simulations. Radiative cooling 
has much stronger effect on the temperature distribution. 
Interestingly, the cooling runs show excess temperature 
at larger radii. We interpret this effect as a consequence 
of the combination of the increased role of shocks in the 
cooling ICM that generate entropy and "adiabatic" com- 
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Fig. 1. — Cross sections through the temperature distribution in the cluster at 2 = 0. All panels correspond to 32h~^ Mpc on a side. 
The top row is for the non-conductive cases while the bottom row is for the runs with anisotropic thermal conduction. Left columns 
corresponds to non-radiative simulations and the right one is for the runs with radiative cooling. The minimum and maximum values of 
the temperature are the same in all panels. This figure illustrates the relative differences between simulations including different physics 
processes. Their effect is quantified in Figure 2 that shows temperature profiles. 



pression fe.g-. lLufkin et al.l(|2000D V Both of these effects 
are more important in the presence of radiative losses in 
the cluster center. The rapid cooling "pulls" the outer 
cluster layers toward the center and heats up the ICM. 
In the presence of cooling the shocks are stronger and 
the heating of the gas more efficient. This effect is illus- 
trated in Figure 3 where we show the entropy profiles. 
The meaning of the curves is the same as in Figure 2. 
Figure 3 shows that the radiative runs exhibit elevated 
entropy at r > 300/i~^kpc. 

We point out that the temperature profile declines with 
the distance from the cluster center even in the runs that 
include conduction. This is possible when the virializa- 



tion shocks heat the gas faster than conduction can re- 
move the heat. 

3.2. Statistical properties of magnetic and velocity field 
orientations 

In the very central parts of the cluster the suppres- 
sion of the effective conduction is further enhanced by 
partially tangential ordering of the fields. This effect is 
shown in the left panel in Figure 4, where we plot the 
absolute value of the radial component of the magnetic 
field vector (left panel) and anisotropy parameter /? for 
the magnetic field (middle panel). This quantity is de- 
fined as /3 = 1 — CTbf/Scr^^, where att and abr are the 
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Fig. 2. — Temperature profiles of the cluster. Dashed green 
curve is for the adiabatic run, solid dark blue for anisotropic 
conduction, dashed light blue for radiative cooling, and solid red 
is for radiative cooling with anisotropic conduction. The profiles 
are not density-weighted. 

transverse and radial magnetic field dispersions, respec- 
tively. For example, abt — (B^)^^^ and the definition 
of abr is analogous. The meaning of the curves is the 
same as in Figure 2. For the absolute value of the radial 
component of the magnetic field, the isotropic case cor- 
responds to 0.5. For magnetic anisotropy, isotropic case 
corresponds to vanishing /3, and tangential and radial 
cases to negative and positive /?, respectively. Horizon- 
tal dotted lines correspond to isotropic configurations. 
Figure 4 demonstrates that the radiative run with field- 
aligned thermal conduction shows a weak bias for tan- 
gential orientation of the magnetic fields. We stress that 
this feature, although technically corresponding to a re- 
gion resolved by up to 16 grid zones, is only tentative 
and that higher resolution runs are required to make any 
definite statements about its nature. In a companion 
paper (Ruszkowski & Oh 2010b, submitted) we make a 
step in this direction by systematically studying the sub- 
structure parameter space and the substructure impact 
on the HBI in isolated cool cores in non-cosmological 
simulations. There is also a tentative decrement in the 
anisotropy parameter at intermediate radii for all four 
runs. This m ay be caused by partial trapping of gravit y 
modes ^Rebusco et al.|[200a [Ruszkowski fc OhI l2010bD . 
At larger radii there appears to be a slight radial bias 
in the orientation of the magnetic field. This could be 
attributed to the accretion along the filaments as the 
magnetic fields are expected t o be locally preferenti ally 
tangential to these structures (jBriiggen et al.|[2C)05b[) or 
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Fig. 3. — Entropy profiles of the cluster. The meaning of 
the curves is the same as in Figure 2. Entropy is defined as 

K = ^'boltzT/n'i^'^ and the plot units are [keV cm^] . 

to inhomogeneous radial fiows in general. The right panel 
of Figure 4 shows the anisotropy parameter /3y for the ve- 
locity field. The definition of this quantity is analogous 
to that for the magnetic fields with velocity dispersions 
replacing the magnetic fields dispersions. However, the 
velocity dispersion is measured with respect to the mean 
streaming velocity of the cluster. This figure bears an in- 
teresting resemblance to its magnetic counterpart. There 
is even stronger tangential bias in the velocity field in the 
radiative run with anisotropic thermal conduction than 
for the magnetic anisotropy parameter for the same run. 
However, there is a significant scatter in these quanti- 
ties and firm conclusion about the trend could be drawn 
from averaging over many clusters. The exact values also 
depend on such factors as the definitions of the cluster 
center and its bulk velocity, and whether the quantities 
are mass- weighted or not. Nevertheless, some similarity 
in these quantities is not unexpected simply due to the 
fact that the magnetic field is frozen to the gas and fol- 
lows it. However, unlike the velocity field, the magnetic 
field has a "memory" of past gas displacement, so the two 
anisotropy parameters are not expected to be identical. 
Intermediate radii tend to have relatively more tangen- 
tial velocity field. 

3.2.1. Radial bias in the velocity and magnetic fields 

As argued i n iSharma et aTl ()2009() and 
[Ruszkowski fc OhI (l2010bD, the ability of the HBI 
and MTI to reconfigure the magnetic field depends on 
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Fig. 4. — Left panel: absolute value of the radial component of the magnetic field vector. Middle and right panels show anisotropy 
parameters for the magnetic (middle panel) and velocity fields. Negative values correspond to tangential fields, zero is for the isotropic 
case. Anisotropy is defined as /3 = 1 — CTj /2(t^ , where rjr and at are the dispersions of either the magnetic field or the velocity field. The 
color coding of all curves is the same as in Figure 2. Horizontal dotted lines correspond to isotropic configurations. 




0.1 1.0 0.1 1.0 

r [Mpc h~^] r [Mpc h~^] 

Fig. 5. — Froude number Fr as a function of the distance from the cluster center (left panel) and the radial component of the gas 
velocity. The radial velocity plot shows net mass-weighted infall velocity. Both panels are plotted using 15 spherical shells of logarithmic 
width 0.168 in the range 15/i~^kpc < r < 5/i~^Mpc. The color coding of the curves is the same as in Figure 2. For Fr > 1 the turbulence 
is expected to begin to randomize the magnetic field configuration that would otherwise be established by the instabilities. 



the level of externally driven turbulence. In the current 
case, the turbulence is driven by the structure formation 
motions. Strong turbulence can significantly affect 
the MTI but some residual r a dial b ias may be present 
in the field. iMcCourt et al.l (|2010t ) showed that even 
for strong turbulence driving, the MTI operates and 
significantly contributes to the velocity power spectrum 
on scales larger than the outer turbulence driving scale 
Lout- On scales smaller than Lout the power spectra in 
the conduction and adiabatic runs are nearly identical. 
As these isolated box experiments considered relatively 
small outer driving scale (Lout = ^pros/16, where /iprcs 
is the pressure scaleheight in a box /ipros/2 high), larger 



outer driving scales may reduce the parameter space 
where the MTI operates. 

If the Froude parameter is less than unity, then the 
turbulence does not reset the magnetic field configu- 
ration established by the MTI. We define the Froude 
number here as Fr = cj/ujbvL, where a is the gas 
velocity dispersion and L is the characteristic length- 
scale. Here we define it entirely for reference purposes 
as the pressure scaleheight L = P/5P in all runs. The 
characteristic frequency uj of the gas perturbed from its 
equilibrium state is 
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dlnS 



d\nr 



dlnT 



d\nr 



(14) 



where g is the total gravity and S = ki,T/n1 is the 
gas entropy. The first term in brackets on the right side 
of Eq. 14 is for the standard Brunt- Vaisala frequency 
while the second term is for its MHD equivalent in 
the presence of anisotropic thermal conduction. The 
Froude number as a function of the distance from the 
cluster center is shown in the left panel in Figure 5. 
The meaning of the curves is the same as in Figure 2. 
For the runs without anisotropic conduction we used 
the standard Brunt- Vaisala frequency to evaluate the 
Froude parameter and for those that include anisotropic 
thermal conduction we used the MHD equivalent of the 
Brunt- Vaisala frequency (see equation 14 above). As 
can be seen in this figure, the Froude parameter typically 
exceeds unity. This suggests that the instabilities should 
be significantly affected by the turbulence. However, 
the estimate of Fr is subject to uncertainties because 
of the arbitrariness of the choice of the characteristic 
turbulence scale. Therefore, Fr parameter should 
only be taken as a rough guide to the impact of the 
turbulence. We also experimented with other definitions 
of L (=100 kpc, hydrostatic pressure scaleheight, and 
iT/|Vaj|, where a is the velocity dispersion and lo is the 
vorticity). In all cases typical values of Fr number tend 
to exceed unity. 

Despite the fact that Froude number tends to ex- 
ceed unity, the velocity and magnetic fields indicate 
clear radial bias for r > 0.8h~^ Mpc in Figure 4. At 
r ~ 0.8h~^ Mpc the result is consistent with no bias 
even though the timescales for the MTI development 
are short (see next section) . Radial bias in the magnetic 
field in the part of the ICM where the temperature 
declines with distance from the center could be caused 
by the MTI. However, the radial bias is present at large 
distances from the cluster center where the contribution 
of turbulence to the pressure support is strongest, which 
tends to isotropize the fields. Thus, a plausible and 
simple explanation for the radial bias is that it is due 
to inhomogeneous radial gas flows either in the bulk of 
the ICM or through the filaments in more distant parts 
of the cluster. We verified that indeed some residual 
net infall velocity is present at larger radii in all four 
runs, i.e., independently of whether anisotropic thermal 
conduction is included. This effect is shown in the right 
panel of Figure 5 where the magnitude of the radial 
component of the gas velocity systematically increases 
at large distances from the center. The meaning of all 
curves is again the same as in Figure 2. This is also 
consistent with the radial bias in the velocity field that 
increases with the distance from the cluster center (right 
panel in Figure 4 discussed above). However, there is 
no excess magnetic and velocity anisotropy bias in the 
runs with anisotropic thermal conduction. Nevertheless, 
we point out that it is difficult to disentangle the 
effect of the MTI and the inhomogeneous radial flow 
on the orientation of the magnetic flelds. Moreover, a 
deflnite statement about the likelihood of the magnetic 
field reorientation due to the MTI in the bulk of the 
cluster volume will have to await a systematic study of 



more than just one cluster as the dynamical states and 
properties of the ICM turbulence (such as the effective 
outer driving scale) will be different in other clusters. 
Here we discuss the result of a pilot study that forms 
the basis for future investigation in this direction. 

The radiative run without conduction (dashed light 
blue line, right panel in Figure 5) exhibits net accretion 
even at smaller radii which is due to the cooling flow. 
It is worth noting that, in the radiative run with 
anisotropic conduction, the weak net infall in the central 
cluster regions is still present but it is significantly re- 
duced compared to the cooling-only case (see solid dark 
blue line for radiative run with conduction and dashed 
light blue line for the radiative run for r < 0.5h~^ Mpc). 



3.2.2. Observations of magnetic field orientation in clusters 

We also note that iPfrommer fc Dursl (|2010| ) reported 
on radio polarization measurements in the Virgo cluster 
based on the observations of magnetic field draping 
around the cluster galaxies. IPfrommer fc Dursl ()2010f ) 
suggests that the polarization pattern in the vicinity 
of the magnetic "tracer" galaxies interacting with the 
ICM is consistent with predominantly radial fields at 
large distances from the cluster center. They further 
pointed out that this observation is consistent with the 
predictions of the MTI. For this mechanism to work, 
the MTI would have to either reestablish the radial 
magnetic field bias on a timescale at least comparable to 
the dynamical timescale of the orbiting "tracer" galaxies 
or the properties of the galaxies must be such that 
they do not dramatically stir the ICM in a continuous 
fashion and/or that a significant time passed since the 
last major merger. Such instability growth rates may be 
possible but whether they are realized depends on the 
details of the dynamical structure of the ICM. 



3.2.3. MTI growth rates and timescales 

In Figure 6 (left panel) we show the MTI growth rates 
in units of Brunt- Vaisala frequency. For the runs with- 
out anisotropic conduction, the rates shown are those 
that would have been present had conduction been in- 
cluded. That is, in this case the MTI rates are com- 
puted via the post processing of the simulation data. In 
all cases, the growth rates were obtained from the MTI 
dispersion relation taking into account the slope of the 
temperature and entropy profiles. As by definition we 
are only interested in those fields that can be reoriented 
in the radial direction due to the MTI, we assumed tan- 
gential fields when computing the MTI growth rates. For 
radial fields the MTI growth rate would be infinite. Fur- 
thermore, we considered perturbation wavelength equal 
to the radius. Shorter fluctuation wavelengths would re- 
sult in somewhat faster growth rates. The right panel in 
Figure 6 shows the corresponding MTI instability growth 
timescales. Both panels in this flgure show that the MTI 
instability has had time to develop across a wide range 
of distances from the cluster center in the non-radiative 
case. Nevertheless, the discussion of the Froude parame- 
ter above suggests that, at least in the cluster considered 
here, the instability might have been overwhelmed by the 
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Fig. 6. — Left: MTI growth rates in units of Brunt- Vaisala frequency. Right: MTI timescales. In both panels the eolor eoding of the 
curves is the same as in Figure 2. For the runs without anisotropic conduction the rates are inferred from the density and temperature 
distributions via the post-processing of the simulation data. 



stirring motions in the ICM. A mild radial bias in the ori- 
entation of magnetic field (and velocity) is present in the 
simulations at large radii but it is consistent with being 
due to either gas inflow or preferentially radial substruc- 
ture motions. However, we reiterate the point made ear- 
lier that it is difficult to disentangle the inhomogeneous 
radial flow from residual MTI. In the runs with radiative 
cooling, the timescales in the more central parts of the 
cluster are long (dashed blue line corresponds to high 
values of the MTI growth timescale for r < 0.2h~^ Mpc 
and the solid red line in the same region is not shown but 
corresponds to very long timescale). This is caused by 
the flatness of the temperature distribution, which makes 
the gas neutrally buoyant. As we explain below, this ef- 
fect leads to an enhanced magnetic field amplification in 
this region when anisotropic conduction is included in 
the simulation. The same stirring motions that tend to 
suppress the MTI in the non-radiative case, are respon- 
sible for driving unimpeded stirring in the radiative one, 
thus amplifying the magnetic fields. However, even in 
the radiative runs, the MTI growth rates are reasonably 
short beyond ~ 0.3h~^ Mpc from the center. We also 
point out that the slopes of the temperature distribution 
in the runs with radiative cooling are steeper at higher 
redshifts, and therefore a wider range of radii was in the 
unstable regime earlier on. 

3.3. Magnetic field amplification 

Cross sections through the cluster center showing the 
distribution of the logarithm of the magnetic field pres- 
sure are shown in Figure 7. The minimum and maxi- 
mum range of magnetic field values is the same in all 
panels. The arrangement of the figures is the same as in 
Figure 1 that shows the temperature distribution: right 
column is for radiative runs and bottom row is for the 



runs with anisotropic thermal conduction. All panels 
show the central region that measures 16 Mpc on the 
side. The magnetic field is clearly amplified in the clus- 
ter center in all four cases. The amplification is much 
stronger in the radiative runs due to the compression of 
the cool gas and magnetic field tha t is f rozen into it. 
This is consistent with the findings of lDubois fc Teyssierl 
(poos) who performed MHD simulations of cluster forma- 
tion with radiative cooling. Interestingly, the radiative 
cooling run with anisotropic thermal conduction (bottom 
right panel) shows even stronger magnetic field amplifi- 
cation than the radiative cooling run. This result is not 
unexpected. The additional amplification occurs in the 
region where the temperature gradient is significantly 
flatter than in other runs. The temperature flattening 
occurs outside the central cool core and up to ~ 0.5h 
Mpc from the cluster center. The nature of convective 
motions changes in the presence of anisotropic thermal 
conduction. When the temperature proflle flattens, the 
ICM tends to become neutrally buoyant. This means 
that the restoring forces in the fluid diminish and the 
substructure infall becomes the main engine for driving 
the unimpeded mixing of the gas. Consequently, this re- 
sults in more efhcient winding up and amplification of 
the magnetic field frozen into the gas via the kinematic 
dynamo effect. 

In order to better quantify the amplification of the 
field, we also plot the distribution of the magnetic field 
along the line passing through the cluster center. This 
is shown in the left panel of Figure 8. The color cod- 
ing of the curves is the same as in Figure 2. The runs 
with cooling boost the field by over two orders of magni- 
tude beyond the amplification seen in the non-radiative 
cases. The top horizontal line denotes the physical field 
at the initial redshift (z — 20) and the bottom one is 
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Fig. 7. — Cross sections through the cluster center showing the distribution of the logarithm of the magnetic field pressure. The 
minimum and maximum range of magnetic field values is the same in all panels. The arrangement of the figures is the same as in Figure 
1 that shows the temperature distribution: right column is for radiative runs and bottom row is for the runs with anisotropic thermal 
conduction. All panels show the central region that measures 8h~^ Mpc on the side. 



for the value of the field that would result from cos- 
mological expansion down to z = without any struc- 
ture formation effects. These reference levels show that 
the magnetic field in clusters is boosted beyond the ini- 
tial physical field by over an order of magnitude. The 
boost in the strength of the magnetic field in the adi- 
abatic case is ^ 10^ compared to the value a uniform 
magnetic field would have at z = 0. The magnitude 
of this ampl i fication is consistent with that obtained by 
IDolag et all (ll999n.lDolag et al.l ()2002D and in our earlier 
work ( Briiggen et al.ll2005bD where the numerical resolu- 
tion was higher. Simple scaling arguments show that 
this field strength exceeds the field expected from the 
magnetic flux freezing arguments. Specifically, for flux 



freezing we have for the ratio of physical field strengths, 
B{z = 0)/Biz = Zini) = [p(z = 0)/pini]2/3. For a typical 
cluster, the final overdensity is ~ 10"^ above the current 
critical density. Thus, for the cosmology relevant to the 
Santa Barbara run we have, B(z — 0)/B{z — Zjni) = 
[10^/(1 -I- Zini)^]^/^ ~ 0.2, which is much smaller than 
the amplification factor of a few x 10 that we see in the 
adiabatic simulation. This shows that the final magnetic 
field in the adiabatic simulation is additionally boosted 
by the shearing motions during the cluster formation. 
We point out that the final magnetic fields at the cluster 
center in the cooling run and the run with anisotropic 
conduction with cooling are boosted even further. The 
magnetic fields are ^ 3 x 10'^ times stronger than the 
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Fig. 8. — The distribution of the magnetic field along the line passing through the cluster center (left panel). The color coding of the 
curves is the same as in Figure 2. The top horizontal line denotes the physical field at the initial redshift (2 = 20) and the bottom one 
is for the value of the field that would result from cosmological expansion down to z = without any structure formation effects. Right 
panel shows the magnetic pressure along the line passing through the cluster. Here the solid red line is for the anisotropic conduction and 
cooling while the dashed light blue line is for the run that includes only radiative cooling. 



physical field at the initial redshift of z = 20 (or equiva- 
lently, over 10^ stronger than a uniform field would have 
at 2 = due to the universal expansion alone, or about 
hundred times stronger than in the adiabatic run) . Such 
levels of magnetic field amplification are also in quali- 
tative agreement with h igher resolution simulations of 
iDubois fc Tevssied ([200l who performed adiabatic and 
radiative runs without thermal conduction. Such field 
strengths are below the MTI and HBI suppression val- 
ues. Thus, for the initial field strengths considered here, 
the suppression of these instabilities by magnetic tension 
does not play a significant role. 

Outside the cluster center the field at z = is reduced 
compared to its initial value. This is simply the result of 
the cosmological expansion. However, the field far from 
the cluster center is still somewhat higher than the uni- 
form field evolved down to 2: = due to the universal 
expansion alone. This could be attributed to mild am- 
plification in filaments or smaller gravitationally bound 
structures. It is interesting to note that the runs with 
radiative cooling show even weaker fields far from the 
cluster center. This is expected as the gas tends to cool 
down and move closer to the center of the cluster, fila- 
ments, and halos amplifying the field in these regions at 
the expense of the rest of the volume. 

As mentioned above, the field amplification in the cool- 
ing run with anisotropic conduction shows stronger field 
amplification compared to its non-conductive counter- 
part (Figure 7, cf. top and bottom panels in the right 
column). In order to better illustrate this effect, we zoom 
in on the central cluster regions and show the magnetic 
pressure along the line passing through the cluster center 
in the right panel of Figure 8. The solid red line is for 
the anisotropic conduction and cooling while the dashed 



light blue line is for the run that includes only radiative 
cooling. It is clear from this figure that the amplifica- 
tion in the former run is stronger by up to two orders of 
magnitude. 

4. CONCLUSIONS 

We performed first magnetohydrodynamical simu- 
lations of cosmological galaxy cluster formation that 
simultaneously include magnetic fields, radiative cooling 
and anisotropic thermal conduction. The presence of 
anisotropic conduction changes the properties of the 
ICM by making it convectively unstable independently 
of the sign of the ICM temperature gradient. In our 
approach, we self-consistently included the amplification 
of the magnetic field due to the shearing motions, 
gas compression enhanced by radiative cooling, and 
the kinematic dynamo associated with the anisotropic 
nature of conduction. Our key findings are: 

1. At large distances turbulent motions tend to reset 
the (radial) orientation of the magnetic field that 
the MTI tries to establish. Nevertheless, some ra- 
dial bias in the orientation of the field is seen at 
large radii in the runs with and without conduc- 
tion. No clear excess of directional bias in the mag- 
netic and velocity field is seen in in the runs with 
anisotropic thermal conduction. The residual bias 
may be due to to the infall of substructures, gas ac- 
cretion along the filaments and the inhomogeneous 
radial fiows in the bulk of the ICM. However, radial 
bias is also expected when the MTI operates. The 
degree of this bias depends on the properties of the 
turbulence driven by structure formation, and dis- 
entangling the inhomogeneous radial gas flows and 
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MTI scenarios is challenging. A systematic the- 
oretical study of a number of clusters is required 
to quantify the role of MTI in shaping the mag- 
netic field topology as the thermal and dynamical 
state of the ICM in other clusters (level of turbulent 
support, MTI growth rates, effective outer scale of 
turbulence driven by structure formation) will vary 
across the mass spectrum of clusters. Here we re- 
port on the first study of a classic Santa Barbara 
cluster that will form the basis for further investi- 
gation to address the above point. 

2. Magnetic field amplification is significantly boosted 
in the presence of radiative cooling which allows 
the gas to concentrate toward the cluster center. 
The central magnetic field at z = is amplified 
by over six orders of magnitude over the value ob- 
tained without substructure formation. 

3. In the presence of anisotropic thermal conduction 
and radiative cooling, the magnetic field is ampli- 
fied by a kinematic dynamo process beyond the 
values obtained in the cooling-only run. This ad- 
ditional amplification occurs over a broad range of 
radii where the temperature profile is relatively fiat. 
In this region, the restoring forces in the fluid di- 
minish and the continuous substructure infall and 
stirring become the main engine for driving the 
unimpeded mixing of the gas and the amplification 
of the magnetic field. 

4. The effective heat conduction from the hotter outer 
layers of the cluster to its center is reduced below 
the full Braginskii-Spitzer value. However, the ef- 
fective radiative cooling driven accretion is notice- 
ably reduced. 

5. The radiative run with anisotropic thermal con- 
duction exhibits a tendency for a tangential bias in 
the velocity and magnetic fields within the inner 
^ 0.3h~^ Mpc that could be associated with HBI 
or trapped gravity modes. Future higher resolu- 
tion cosmological cluster formation studies will 
assess the robustness of this effect. Any possible 
anisotropy in the magnetic field distribution may 
be detectable via radio polarization measurements 
with Square Kilometer Array and LOFAR, while 
the bias in the velocity field could be probed 
with the future International X-ray Observatory 
mission. 
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5. APPENDIX A 

In order to check the implementation of the cosmo- 
logical terms in the MHD equations we evolved spatially 
constant matter density and magnetic field while neglect- 
ing any velocity perturbations. The result of this test is 
shown in Figure 9. In the left panel we show the scaling 
of the temperature (in arbitrary units) with the cosmo- 
logical expansion factor a. Shown are the code result 
for the physical temperature and the power law fit (solid 
and dashed lines are practically indistinguishable). The 
slope of T{a) oc agrees with the standard theoretical 
expectation. The right panel shows the physical mag- 
netic field in arbitrary units as a function of the scaling 
parameter a and a powerlaw fit to this relation, both as a 
function of the scaling parameter a. Here again the fit is 
perfect, the solid and dashed lines overlap, and the slope 
of the field is {B{a) oc a^^) ensures the conservation of 
the magnetic fiux. 

6. APPENDIX B 

In order to test the implementation of the anisotropic 
thermal conduction module we compared linear theory 
MTI growth rates with the code results. This test is very 
similar to the one discussed in Parrish & Stone (2008). 
That is, we set up a two dimensional stratified hydro- 
static atmosphere with very shallow density and temper- 
ature profiles such that 

Tiz) = To{l-y/yo) (15) 



p(z)=po(l-y/y„)2, (16) 

where To and po are constants. The characteristic length- 
scale yo was set to 1% of the horizontal height of the 
computational box. The gravitational field was assumed 
constant throughout the computational domain. We set 
hydrostatic boundary conditions on the top and bottom 
boundary. We also impose constant temperature in the 
boundary zones to prevent the escape of thermal energy 
from the simulation box via thermal conduction. That 
is, the computational domain is effectively adiabatic. We 
set periodic boundary conditions in the horizontal direc- 
tion. Initial magnetic field is horizontal with very high 
magnetic f3 parameter. In order to seed the instability we 
introduce a very small sinusoidal velocity perturbation 
in the vertical direction. Specifically, the initial vertical 
component of the velocity field is given by 

v,{y) = VoSm{2TTy/L), (17) 

where Vo is a very small constant amplitude of the veloc- 
ity perturbation and L is the box size in the horizontal 
j/-direction. We evolve such initial conditions and, using 
the velocity field as a function of time, we compute the in- 
stabilit y growth rate a fo l lowin g a prescription similar to 
that in iParrish fc Stoiii (|2005f ). These results are then 
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Fig. 9. — The evolution of physical tempcratairc (left panel) and physical magnetic field with the cosmological expansion parameter a. 
Quantities on vertical axes in both panels are in arbitrary units. The test assumes uniform magnetic field, density and temperature and 
neglects any seed fluctuations in the velocity field. The results are in perfect agreement wit the theoretical scaling relations (powerlaw fits 
to the computed relations are also plotted and match these relations perfectly.) 



compared to the linear theory prediction in the nondi- 
mentional form 

+ -a^x + (T+ -jT—7;X = 0, (18) 
7 a In 6 



where k is the wavenumber of the velocity fluctuation. 
The result of this test is shown in Figure 10. The growth 
rates predicted by the code are in excellent agreement 
with the theoretical expectations. 



where cr is the instability growth rate in units of the 
Brunt- Vaisala frequency ojbv and 



x = h- l)K-euJsv^ (19) 
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Fig. 10. Magnctothcrmal instability growth rate 
in units of Brunt- Vaisala frequency as a function of 
a; = (7 — 1)k{T / P)k'^ /lobv t where k is the wavenumbcr of 
the velocity fluctuation and lobv is the Brunt- Vaisala frequency. 
Solid line denotes the prediction from linear theory and the dots 
show the results obtained with the FLASH code. 
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